Implicit-explicit (IMEX) evolution of single black holes 
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Numerical simulations of binary black holes — an important predictive tool for the detection of 
gravitational waves — are computationally expensive, especially for binaries with high mass ratios 
or with rapidly spinning constituent holes. Existing codes for evolving binary black holes rely on 
explicit timestepping methods, for which the timestep size is limited by the smallest spatial scale 
through the Courant-Friedrichs-Lewy condition. Binary inspiral typically involves spatial scales 
(the spatial resolution required by a small or rapidly spinning hole) which are orders of magnitude 
smaller than the relevant (orbital, precession, and radiation-reaction) timescales characterizing the 
inspiral. Therefore, in explicit evolutions of binary black holes, the timestep size is typically orders 
of magnitude smaller than the relevant physical timescales. Implicit timestepping methods allow 
for larger timesteps, and they often reduce the total computational cost (without significant loss 
of accuracy) for problems dominated by spatial rather than temporal error, such as for binary- 
black-hole inspiral in corotating coordinates. However, fully implicit methods can be difficult to 
implement for nonlinear evolution systems like the Einstein equations. Therefore, in this paper 
we explore implicit-explicit (IMEX) methods and use them for the first time to evolve black-hole 
spacetimes. Specifically, as a first step toward IMEX evolution of a full binary-black-hole spacetime, 
we develop an IMEX algorithm for the generalized harmonic formulation of the Einstein equations 
and use this algorithm to evolve stationary and perturbed single-black-hole spacetimes. Numerical 
experiments explore the stability and computational efficiency of our method. 
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I. INTRODUCTION 

Binary black holes (BBHs) are important sources of 
gravitational waves for the current and future gravita- 
tional wave detectors such as LIGO, Virgo, LCGT [lHl] 
and LISA Data-analysis of these gravitational wave 

detectors proceeds with matched filtering, which requires 
accurate knowledge of the expected waveforms. This mo- 
tivates numerical simulations of the inspiral, merger and 
ringdown of two black holes. Starting with Pretorius' 
2005 breakthrough 0], several research groups have de- 
veloped numerical codes capable of simulating this pro- 
cess (see [|[ for a recent review) . 

BBH inspiral simulations for gravitational wave detec- 
tors must cover at least the last « 10 orbits of the inspi- 
ral, and possibly many more [ij-Hsj , requiring simulations 
significantly longer than the dynamical timescales of the 
individual black holes. This separation of temporal scales 
becomes particularly pronounced for a BBH with mass- 
ratio g > 1: The dynamical time of the smaller black 
hole shrinks proportional to 1/q. Simultaneously, the in- 
spiral proceeds slower and the time the binary spends in 
the strong-field regime lengthens proportionally to q. 

All published numerical simulations of BBH inspi- 
ral and merger employ explicit timestepping algorithms 
which are subject to the Courant-Friedrichs-Lewy (CFL) 
condition which limits the timestep size by the small- 
est spatial scale in the problem. Binary inspiral typ- 
ically involves spatial scales (the spatial resolution re- 
quired by a small or rapidly spinning hole) which are 



orders of magnitude smaller than the relevant (orbital, 
precession, and radiation-reaction) timescales character- 
izing the inspiral. In explicit binary evolutions the CFL 
condition then effectively fixes the timestep size to be the 
dynamical timescale (see the last paragraph) for one of 
the constituent holes. Such a timestep is orders of mag- 
nitude smaller than the relevant physical timescales for 
the binary as a whole; particularly when the binary has a 
large mass ratio (such as the simulations in Refs. (l^.[l5|) 
or when at least one constituent hole has a high spin 
(since the horizon of the high-spin hole then requires 
higher spatial resolution). For instance, a simulation 
with constituent holes with dimensionless spin magni- 
tudes 0.95 [l(| required half a million timesteps over 12.5 
orbits. 

Were the CFL restriction overcome, computation of 
BBH inspirals with higher mass ratios, higher spins, and 
more orbits could become feasible. Implicit timestep- 
ping is one way to overcome the CFL condition and 
take larger timesteps. Of course, larger timesteps cor- 
respond to larger temporal truncation errors; however, 
a small timestep is required in BBH inspirals for stabil- 
ity (CFL condition) rather than accuracy (since, as ar- 
gued above, the accuracy of a BBH inspiral is typically 
limited by spatial resolution, not temporal resolution). 
For problems dominated by spatial rather then tempo- 
ral error, implicit timestepping methods often reduce the 
total computational cost (without significant loss of ac- 
curacy), but fully implicit methods can be difficult to im- 
plement for nonlinear evolution systems like the Einstein 
equations. Implicit-explicit (IMEX) methods [l7l - |2"oj are 
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a compromise which we explore here. IMEX timestep- 
ping has been successfully applied to a variety of prob- 
lems, including fluid-structure interaction [Hj], relativis- 
tic plasma astrophysics {22j . and hydrodynamics with 
heat conduction 23 1. In Ref. J24[, Lau, Pfeiffer, and Hes- 
thaven applied IMEX methods to evolve a forced scalar 
wave propagating on a curved spacetime (a Schwarzschild 
black hole), achieving stable evolutions with timcstcp 
sizes « 1000 times larger than with explicit methods. 

In this paper, we lay much of the groundwork toward 
applying IMEX methods to full binary-black-hole evolu- 
tions. We develop an IMEX algorithm for one partic- 
ular formulation of Einstein's equations used in explicit 
BBH evolutions, the generalized harmonic formulation 
(see [25| and references therein). We use our IMEX 
algorithm to perform the first IMEX evolutions of sin- 
gle black holes (both static and dynamically perturbed) . 
Our single-black-hole evolutions demonstrate the stabil- 
ity of our IMEX method. Further numerical experiments 
also investigate our method's efficiency; the IMEX algo- 
rithm offers a computational cost competitive with ex- 
plicit evolution for sufficiently large step sizes. (Note that 
improved efficiency docs not automatically follow from an 
IMEX algorithm affording larger timcsteps, since each 
IMEX timestep is more expensive than an explicit step.) 
We also discuss further efficiency improvements of our 
IMEX implementation, and provide an outlook toward 
simulation of black hole binaries with IMEX techniques. 

This paper is organized as follows. In Sec. [HI we derive 
the IMEX generalized harmonic equations and boundary 
conditions that we will use. In Sec. lIIII we explore numer- 
ical simulations using these equations, with a particular 
focus on the stability and efficiency gains of these simu- 
lations. We conclude in Sec. II VI by discussing the impli- 
cations of our results, emphasizing the probable gains in 
computational efficiency when using IMEX in full binary- 
black-holc simulations. 



II. IMEX FORMULATION OF EINSTEIN'S 
EQUATIONS 

The generalized harmonic formulation of Einstein's 
equations consists of ten coupled scalar wave equations. 
Therefore, the present discussion will borrow heavily 
from our earlier work on IMEX evolutions of scalar fields 
on curved backgrounds (24|. 



of Ref. [25[) is the following: 

dt^ab = (1 + li)V k d k iP ab - NTl ab - jiVHkat, (la) 

d t u ab = v k d k n ab - Ng jk d^ kab + w 2 v k d k ip ab 

+ 2Nii cd (gi k $ jca $kd b - n cQ rU - i> ef T ace T bd} ) 

- 2NV [a H b) - \Nt c t d Il cd Il ab - Nt c U CJ g lk <$> kab 

+ l0 N(2S c {a t b) ~ ^ ab t c ) (H c + r c ) - lll2 V k <5> kab 

(lb) 

dt^ jab = v k d k $ jab - Nd,n ab + N l2 d 3 tp ab 

+ lNt c t d $ jcd Il ab + Ng km t c $ jkc $ mab - N~f2*jab. 

(lc) 

Here, N, V k , and gj k are the spacetime metric's associ- 
ated lapse function, shift vector, and spatial metric in- 
duced on level-i slices. Latin indices from the middle of 
the alphabet = 1,2,3 range only over spatial di- 

mensions. As a one-form, t a = —Nd a t is the unit normal 
to the temporal foliation defined by the coordinate time 
t. The other fundamental variables II a f, = —t c d c ipab and 
Q kab = d k ipab arise from the reduction of the general- 
ized harmonic equations to first order form. The latter 
definition leads to the auxiliary constraint 



C ka b = d k 1pab - ®kab = 0. 



(2) 



The variable r a = ip bc T abc represents a contraction of 
the Christoffel symbols T abc of the spacetime metric tp ab . 
Time derivatives dtipab inside T abc are evaluated in terms 
of N, V k ,U ab , and $ fcab Hg. 

The functions H c are freely specifiable and embody 
the coordinate- freedom of Einstein's equations [25J . Ein- 
stein's equations can be written as a set of constrained 
evolution equations; in the generalized harmonic formu- 
lation, the fundamental constraint takes the form 



C a = H a + F 



0. 



(3) 



Constraint damping 0, [25l - |27| is used to enforce both 
the fundamental constraint ([3]) and the auxiliary con- 
straint ([2]). Those terms in Eqs. ([T]) proportional to 

70 damp the fundamental constraint @. Those terms 
proportional to 71 and 72 in Eqs. (JTJ) damp the con- 
straint (|2J|. Our IMEX formulation converts to second 
order variables and so the auxiliary constraint is triv- 
ially satisfied. Therefore, in the rest of this paper, we set 

71 = = 72 in all IMEX evolutions. 



A. Generalized harmonic system 

Our goal is to solve Einstein's equations for the space- 
time metric ip a b, where Latin indices from the start of the 
alphabet (a, . . . , /) range over 0, 1, 2, 3. The first order 
generalized harmonic formulation of the Einstein evolu- 
tion equations given by Lindblom et al (Eqs. (35)-(37) 



B. First-order implicit equations and second-order 
implicit equation for the metric 

Although ([1]) is a system of partial differential equa- 
tions (PDEs), we formally view it as an ordinary differ- 
ential equation (ODE) initial value problem, 

du 

— = f(t,u), u{t )=u Q , (4) 
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so that our notation conforms with the literature 
[20[ on IMEX ODE methods. [Otherwise, we would have 
used partial time differentiation in (J4j) .] The system ([T]) is 
actually also solved as an initial boundary value problem; 
however, we defer the issue of boundary conditions to a 
later subsection. In this view u represents the collection 
(ipab, n a b, $fe a b) of fundamental fields. Furthermore, we 
assume there exists a splitting 



f{t 1 u) = f(t,u) + f E {t,u) 



(5) 



and 



of the right-hand side / into an explicit sector 
an implicit sector f 1 . In this paper, as in Ref. [24(, we 
split by equation. That is, we choose which terms on the 
right-hand side of Eq. ([T]) are to be treated implicitly. 

To take a timestep, we choose an IMEX timestepping 
algorithm, such as ImexEuler, Additive Runge Kutta 
(ARK) [23, or semi-implicit spectral-deferred correction 
(SISDC) |17M20j . We note that while ARK was used 
almost exclusively in Ref. [24[, we have encountered sta- 
bility issues with its use in the work presented here, and 
therefore focus here on SISDC. As explained in Sec. II A 
of Ref. , each of these algorithms requires that we are 
able to solve (multiple times per timestep) an implicit 
equation of the form 



u — af 1 (t, u) = B, 



(6) 



where a is proportional to the step size At and the inho- 
mogeneity B is defined by the algorithm. For example, 
the corresponding equation for ImexEuler integration, 

u n+1 - Mf(t n+1 ,u n+ i) = u„ + Atf E (t n ,u n ), (7) 

is solved to advance the solution from time t n to time 
t n+ i. Concrete expressions for B are given in Ref. [24| 
for ARK and in Appendix U for SISDC. 

The IMEX splitting of the system ([1} that we chose is 
analogous to the "case (ii)" equations for the scalar- wave 
system given as Eqs. (15a)-(15c) in Ref. [Hj]. Specifi- 
cally, we treat implicitly the entire right-hand sides of 
Eqs. (JTJi,) and (JlJ;). However, a fully implicit treatment 
of the equation for II a h has turned out to be prohibitively 
complicated. Therefore, of the terms appearing in the 
right-hand side of Eq. ([lb), we have chosen to include 
in the implicit sector only the principal-part terms and, 
possibly, the constraint damping term proportional to 70 . 
The principal-part terms are the stiff terms which most 
constrain the timestep size, and, as we shall see later, the 
constraint damping term is also stiff. Implicit treatment 
of the remaining terms on the right-hand side of Eq. (JTJd) 
would be difficult because the implicit equation which re- 
sults from their inclusion has an extremely complicated 
variation. This variation would be required were the re- 
sulting equation solved (as part of the overall system) via 
Newton iteration. 

Our splitting of Eq. ([T}d) could be improved upon. In- 
deed, with fu ab (t, u) representing the right-hand side 
of the evolution equation (TT]d) for Tl ab , a binary evolu- 
tion based on the dual-frames approach will have fn ab = 



O(w), where to is the orbital frequency (a small quantity). 
However, for our described splitting both /£ and 
would be O(l). Although their combination is small, 
each individual term on the right-hand side of (QJi) need 
not be. In other words, there appears to be no natu- 
ral splitting by equation for Eqs. ([T]), as there often is 
for, say, advection-diffusion problems. While we do not 
yet fully appreciate the consequences of the splitting we 
shall employ here, we are considering approaches to mit- 
igate potential problems with our splitting-by-equation 
approach. Among these is a fully implicit implementa- 
tion of Eq. ((T}d) , with other possibilities discussed in the 
conclusion of Ref. [24[ ■ 

Our choices above correspond to the following first- 
order implicit equation for H a b'- 



TL ab - a[V k d k U ab - Ng jk dj^ kab 

+ Y N(26 c {a t b) - t/W c ) (H c + r c )] 



B, 



(8) 



Here we have split the damping parameter as 70 = 7o + 
7o\ which in general allows for part of the damping term 
to be treated implicitly (if 7q ^ 0) and part explicitly (if 
7(f ^ 0). In Eq. JgJ) we view T e as 



r, = 



dn cd 



n 



cd 



r e — — —lied 



(9) 



terms with n a6 terms without U a 



with the details of this decomposition given in Ap- 
pendix [A] The reason for the decomposition is given 
immediately after Eq. (fTB")) . In all, our first-order im- 
plicit equations corresponding to the evolution system 
([T]) are then as follows: 

^ ab - a(V k d k i> ab - NU ab ) = B,^ ah (10a) 
n Qfa - a (v k d k n ab - A '/'^ <!>■„, + NQ ab cd U cd 

+ NG ab ) = B Uab (10b) 

$ jab - 'I (I ' 'I',,,;, - NdjIU + \Nt c t d ^ 3cd Ii ab 



Ng km t c <S> ]kc <S> 



jnab 



= B« 



(10c) 



where 



Qab cd = ll{28\ a t b) ~^ ab t e 



£>£e 

an 



g ab = y (25\ a t b) -i> ab t e ) 



cd 



-n 



cd 



(ii) 



To solve these equations, we first take a combination 
of them to get a single second-order equation for ip ab . In 
terms of £ ab = ip ab — aV k dkip a b, we express (flUk ) as 



aNIl 



a!) 



B,i 



£,ab 



(12) 



Multiplication of Eq. (|10b ) by aN, followed by a substi- 
tution with ((T2|) , yields 



aATI afc - a 2 NV h d k Il ab + a 2 N 2 g^d^ kab 
- aNQ ab cd {B. 4 , cd - - a 2 N 2 g ab = aNBn ab . (13) 
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The decomposition J9]) ensures that the substitution with 
Eq. (|T2"j) is also made for the H cd terms in T e . We subtract 
the last equation from (fTOk ) to reach 

Vab - aV k d k ip ab + a 2 NV k d k Tl ab - a 2 N 2 g jk d^ kab 
- aNQ ab c % d + a 2 N 2 G ab 

= B i , ab -aNB Ilab -aNQ ab cd B iJcd . (14) 

We must eliminate the term a 2 NV k d k H ab from the re- 
sult. To this end, we contract Eq. (fTOb ) into aV j , thereby 
finding 

aV j $ jab - a 2 V k V j d k $ jab + a 2 NV J d 3 IV ab 

- ±a 2 Nt c t d V j $ jcd Il ab - a 2 Ng km t c V j $ jkc $ mab 

= aV j B 9jab , (15) 

which, using Eq. (fl~2)) . we rewrite as 

aV j $ jab - a 2 V k V'd k ^ ab + aPNVtdjlleb 
+ \at c t d V^ ]cd U - a 2 Ng km t c V j $ jkc $ mab 
= aV J B 9jab + lat c t d V j ^ jcd B^ ab . (16) 

Subtracting the last equation from (|14[) and making sub- 
stitutions with the constraint ([2]), we arrive at the fol- 
lowing second-order equation: 

4> ab - 2aV k d k ^ ab - a 2 (N 2 gi k - \ ' V 1 )il,ih >:,>, 

- \at c t d V={d^ cd )^ ab ~ aV k d k 4, ab ) 
+ a 2 Ng km t c V 1 {d^ kc ){d m il, ab ) 

- aNQ ab cd {^j cd - aV k d k ^ cd ) + a 2 N 2 Q ab 

= (1 - \at c t d V'd^ cd )B^ b - aNBn ab - aV k B^, kab 

- aNQ ab cd Bjp d + terms homogeneous in C ka b- 

(17) 

To solve the system (|T0|) . we first solve (fT7|) . subject 
to boundary conditions discussed in Sec. Ill CI Next, 
we recover H ab algebraically from (fTUk). Finally, we 
set $fc a fc = d k ip a b, i-e., we enforce that the constraint 
C ka b = 0. 

We stress that, as a linear and undifferentiated com- 
bination of Eqs. (fTU| for the first-order system, Eq. p7|) 
actually contains no second-order derivatives of ip ab . In- 
deed, all of the i3-terms on the right-hand side of Eq. p7|) 
appear undifferentiated, indicating that we have not dif- 
ferentiated the first-order system (|10[) . Each second- 
order derivative of t/) ab on the left-hand side of (jl"7|) is 
precisely canceled by a corresponding term appearing 
in one of the constraint terms on the right-hand side 
[not shown explicitly in Eq. (|17[) ]. Now, when numer- 
ically solving Eq. (|1T[) . we set the constraint terms from 
the right-hand side to zero, thereby creating a genuinely 
second-order equation. We discuss the permissibility of 
this procedure in Sec. Ill Dl below. 



C. Boundary conditions 

For black-hole evolutions which employ excision, the 
inner boundary lies within an apparent horizon. For this 
scenario we adopt no inner boundary condition, regard- 
less of what condition is adopted at the outer boundary 
and despite the fact that Eq. (fTTf is a second-order equa- 
tion. In the context of scalar fields on a fixed black-hole 
background, Ref. [24| has discussed the motivation for 
and permissibility of this procedure. A similar analytical 
treatment of the coupled nonlinear system (|17|) would be, 
we suspect, a difficult piece of mathematical analysis, one 
beyond the scope of this paper. Therefore, here we con- 
tent ourselves both with the scalar field analogy and the 
observation that the lack of an inner boundary condition 
has caused no difficulties numerically. Nevertheless, the 
issue merits further study. 

The outer boundary condition that we apply to 
Eq. (fT7)) is either (i) a fixed Dirichlct condition on each 
component ip ab of the spacetime metric or (ii) the fol- 
lowing condition. In terms of the incoming characteris- 
tic variable U~ b = H ab — n k $ kab (where n k is the unit, 
outward-pointing, normal vector to the boundary), we 
rewrite Eq. (fTUk ) as 

i' ab + a(Nn k - V k )d k ip ab = B^ ab - aNU~ b + aNn k C kabl 

(18) 

We control U ab at the boundary; therefore, both B^ ab 
and U~ b here appear as fixed quantities, and Eq. (|T8]) 
represents a boundary condition on ip ab . Moreover, when 
numerically enforcing this condition we also set the con- 
straint term on the right-hand side to zero. 



D. Implicit equation for the auxiliary constraint 

Eqs. (fTUa ) and (fTUb ) imply an implicit equation for 
the auxiliary constraint. Partial differentiation of (fTOk ) 
yields 

d^ ab - a[(d V k ){d k i> ab ) + v k d k d^ ab 

- (djN)IL ab - NdjILab) = djB^ ab . (19) 

To express the derivatives of the lapse and shift in terms 
of derivatives of the metric ipab, we use the result 

6<p ab = -2N-H a t b 6N - 2N~ 1 g k{a t b) 5V k + g\ a g k b) 5g lkl 

(20) 

which in turn yields 

6N = ~±Nt c t d 6^ cd , SV k = Ng km t c H mc . (21) 

Insertion of these results (with the variation 5 — > dj) into 
(H3) gives 

d^ ab -a[Ng km t c (d 3 i; mc )(d k ^ ab ) + V k 8 k d^ ab 

+ \Nt c t d {d^ cd )U ab - NdjUab] = d 3 Bip ab . (22) 
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Finally, we subtract (fTUb ) from the last equation and 
make substitutions with the constraint to reach 

Cjab - a[V k d k C jab + Ng km t c (<S> ]kc C ma b+C jkc d m i>ab) 
+ \Nt c t d Cj C dIi a b\ = djB^ ab — B^ jab . (23) 

This equation is analogous to Eq. (20) of Ref. [2~H . 

C 3 - a£ v Cj = djB^, - B*, , (24) 

for scalar waves on a fixed curved background, where 
the overbar on Cj serves to differentiate this constraint 
from the generalized harmonic constraint C a in Eq. ([3]) 
(which carries a spacetime rather than spatial index in 
any case). Specifically, in the scalar wave scenario the 
variables (ip, II, $fe) are analogous to the generalized har- 
monic variables (ij) a b, n o j,, ^kab)j an d the auxiliary con- 
straint is Cj = djtp — Starting with a prescribed Cj at 
the outer boundary, we may integrate Eq. ([24]) along the 
integral curves of the shift vector. This independent inte- 
gration of Cj proved important toward understanding in 
what sense solving the second-order implicit equation for 
tfj [analogous to Eq. (fTT)) ] was equivalent to solving the 
first-order system for (-0,11, $fe) [analogous to Eq. (fpQ|) ]. 
Such an independent integration of (|23[) is clearly not 
possible. Nevertheless, provided both Cj a t = on the 
outer boundary and a vanishing right-hand source in 
(|23p . the equation formally determines Cj a b = along 
the integral curves of V k . This motivates our neglecting 
the terms homogeneous in Cj a b in Eq. (fT7|) . 

Consideration of our steps above for solving (flU)) shows 
that the constraint Cj a b remains exactly zero throughout 
our IMEX scheme. We are then effectively evolving only 
the variables ip a b and H a b- Our reasons for nevertheless 
retaining &j a b in the formalism are twofold. First, SpEC 
— the software project we have used for simulations — 
chiefly supports first order symmetric hyperbolic sys- 
tems. Second, as described in the conclusion, for the 
binary problem we envision a split by region approach, 
in which outer subdomains are treated explicitly and in- 
ner subdomains (spherical shells) immediately near the 
holes are treated by IMEX methods. Since explicit evo- 
lutions in SpEC currently require a first order system, the 
variable §j a b must be present in the outer subdomains. 
Coupling between the outer and inner subdomains is then 
facilitated by having ( I> ; „; ; also available on the inner sub- 
domains. There has been recent progress in applying 
spectral methods to evolve second order in space partial 
differential equations [29j|. If these techniques work for 
the generalized harmonic system, it should be possible to 
abandon &j a b entirely. 

III. NUMERICAL EXPERIMENTS 

Through numerical simulations of single black holes, 
we now examine the behavior of the scheme presented 
above. We evolve initial data representing both (i) the 
static Schwarzschild solution in Kcrr-Schild coordinates 



and (ii) the same solution with a superposed ingoing 
pulse of gravitational radiation. The latter is a vac- 
uum problem with non-trivial evolution. As the gravi- 
tational wave pulse travels inward, it hits and perturbs 
the black hole. Most of the pulse is absorbed by the black 
hole, increasing its mass; the rest is scattered and prop- 
agates away. This test features initial dynamics on short 
timcscales (moving pulse of radiation, perturbed black 
hole), with relaxation to time-independence. Eventually, 
the black hole settles down to a stationary black hole, 
and the scattered radiation leaves the computational do- 
main through the outer boundary. Technical details for 
the dynamical case (ii) are summarized in Appendix [Cj 



A. Long-time stability of IMEX evolutions 

In this subsection we demonstrate the stability of our 
IMEX algorithm by evolving the static Schwarzschild so- 
lution in Kerr-Schild coordinates to late times (up to 
10 4 M), adopting fixed Dirichlet conditions, that is with 
tjj a b fixed as the analytical solution on the outer bound- 
ary. We note that the radiation conditions (fT5|) . with 
U~ b determined by the analytical solution on the outer 
boundary, apparently give rise to an extremely weak in- 
stability. Indeed, with Eq. (fT5)) a slowly growing insta- 
bility appears after (sometimes well after) time 10 3 M. 
We specify no inner boundary condition (cf. Sec. Ill Cp . 
Our domain, a single spherical shell with Cartesian cen- 
ter (0.01,-0.0097,0.003), is determined by a top spher- 
ical harmonic index ^ max — 7 and the radial interval 
1.9 < r < 11.9, with N r = 15 radial collocation points 
and an exponential mapping of the radial coordinate (see 
Eq. (48) of HH). Results for Cartesian center (0, 0, 0) are 
qualitatively similar, but with the corresponding errors a 
few orders of magnitude smaller. For constraint damping 
parameters, we have taken 7q = 1 and Jq = 0. 

We have performed IMEX evolutions with an 
ImcxEulcr timestepper (first order accurate and requir- 
ing one solution of the system (JTUJ) per timestep), 3-point 
(substep) Gauss-Lobatto SISDC (GLoSISDC3, fourth or- 
der accurate, eight implicit solves per timestep), and 2- 
point (substep) Gauss-Radau-right SISDC (GRrSISDC2, 
third order accurate, six implicit solves per timestep). 
Since the geometry is time-independent, numerical solu- 
tion of (|17[) will be achieved without any iterations in the 
Newton-Raphson algorithm, assuming that the solution 
at the previous timestep serves as an initial guess. To pre- 
vent this trivial convergence, we have rcscalcd the initial 
guess ip® b — > 1. 00001 V^;, before each implicit solve. For 
GLoSISDC3 and GRrSISDC2 respectively, Figs. □ and 
[5] depict error histories for the metric ip a b as measured 
against the exact solution. Each plot exhibits long-time 
stability for the larger timesteps considered but weak in- 
stability for some of the smaller timesteps. 

Examination of the stability diagrams for these meth- 
ods suggests a heuristic explanation of our results. The 
diagram for a given (either explicit or implicit) ODE 
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FIG. 1. Error histories for GLoSISDC3. || • || represents the 1- 
norm with respect to the Cartesian coordinate measure over 
the spherical shell, i. e. ||/|| = f v \f\dxdydz, and V is the 
improper (coordinate) volume of the spherical shell. As men- 
tioned in the text, Aip a b denotes the difference between the 
numerical metric and the exact solution. 
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FIG. 3. Diagram for implicit sector of GLoSISDC3. The 
bottom plot depicts the cross section of the top plot along 
the imaginary axis, with A = 177 6 iR. 
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FIG. 2. Error histories for GRrSISDC2. See the caption of 
Figure [TJ for an explanation of the figure labels. 



method is determined by its application to the model 
problem du/dt = Xu, where A = £ + vq. Subject to the 
initial condition Uq = 1, a single timcstcp for a given 
method produces an update UAt = Amp(AAi), the am- 
plification factor which is a function of the complex vari- 
able XAt. The region of absolute stability for a given 
method is then the domain in the (AAi)-plane for which 
|Amp(AAf)| < 1. Figures [3] and 0] respectively depict 
the stability diagrams for GLoSISDC3 and GRrSISDC2, 
with the model problem treated fully implicitly, i.e. with 
f 1 = Xu and f E = 0. For both diagrams, our interest lies 
with the imaginary axis, since the system ([1} of equations 
we evolve supports the propagation of waves. 

For GLoSISDC3, the imaginary axis lies within the re- 



gion of absolute stability, except for a portion around 
the origin. The bottom panel of Fig. |3] shows that 
|Amp(iryAf)| > 1 for \r)At\ < 1.28, with the maximum at 
rjAt w ±1. Note also that \r]At\ < 0.35 corresponds to an 
essentially conservative method, since then |Amp(i?7At)| 
is very close to unity. Therefore, assuming A in the model 
problem is purely imaginary, we expect growth in the 
numerical solution for timesteps At < 1 . 28 1 A | 1 , and ab- 
solute stability for At > 1 . 28 1 A | 1 . Figure [4] provides 
the analogous information for GRrSISDC2; the bottom 
plot indicates growth for timesteps At < 0.51 1 A| ~ x but 
absolute stability for At > 0.51|A| _1 . We now attempt 
to identify A in the model problem with characteristic 
speeds for the evolution system (p}. 

Given an outward-pointing unit normal n (often to 
the boundary of a computational domain or subdomain) , 
the characteristic variables of Eqs. ([1]) are 



ipab, n a6 ±n fe $ 



kab 1 



(<J| - njn k )$ kab , 



(25) 



and their respective characteristic speeds are 



n k V k , -n k V k ±N, -n k V k . (26) 



Equations (J23J> and are derived in [25( [see Eqs. (32)- 
(34) of that reference and the text thereafter, but set 
72 = = 71 as is the case here]. For the Schwarzschild 
solution in Kerr-Schild coordinates (see Eq. (34) of [Hj]), 
the characteristic speeds for propagation orthogonal to 
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FIG. 4. Diagram for implicit sector of GRrSISDC2. See rele- 
vant comments given in the caption of Fig. [3] 



an r = const sphere reduce to 



n k V k 



n k V k ±N = 



2M 



\A" 2 + 2 Mr 

2M 
Vr 2 + 2Mr 



± 



r + 2M 



(27a) 
(27b) 



where these expressions correspond to coordinate spheres 
adapted to the spherical symmetry, i.e. to Cartesian cen- 
ter (0,0,0). The smallest speeds (in magnitude) are 
n] i V k near the outer boundary (r large), and n k V k — N 
near the horizon (r = 2M). 

An instability driven by the speed Eq. (|27a[) evaluated 
at the outer boundary appears consistent with the stabil- 
ity diagrams Figs. [3] and @] in the following sense: At the 
outer boundary r = 11.9, n k V ~ 0.16. Assuming wave 
solutions propagating with this characteristic speed, we 
have A = i0.16 in the model problem above. Our sim- 
ple analysis predicts instability when At < 8.0 for GLo- 
SISDC3 and At < 3.2 for GRrSISDC2, with stability for 
At larger than these estimates. The results depicted in 
Figs. [T] and [2] arc consistent with these predictions. 

Note that the bottom panels of Figs. [3] and @] indi- 
cate better stability properties for \r)At\ close to zero. 
However, even if the characteristic speeds at the outer 
boundary correspond to this "near-stable" portion of 
the imaginary axis in the relevant stability diagram, the 
characteristic speeds normal to r = const surfaces for 
smaller radius r have larger characteristic speeds, and 
thus |Amp(i7yAt)| near its maximum. 

Moreover, the predictions of our stability analysis ap- 
pear at least qualitatively correct when the location of the 
outer boundary is moved to larger radii, where n k V k is 



smaller. As rikV k decreases, larger timcsteps At should 
become unstable. Indeed, with GLoSISDC3 for exam- 
ple, we find that At = 8 is unstable for r = 18.9 (and 
apparently independent of radial resolution). By simi- 
larly pushing the outer boundary outward, we can render 
At = 4 unstable for GRrSISDC2. Finally, we note that 
the standard stability region for backward Euler contains 
the entire imaginary axis, and is dissipative for imaginary 
A. All of our evolutions with ImcxEuler have proved 
correspondingly stable, even for small timesteps (with 
At — 1/2 the smallest considered). 



B. Convergence of the IMEX method 

We now verify both the temporal and spatial conver- 
gence of our scheme, using the perturbed initial data 
[case (ii)] described both above and in more detail in 
Appendix O We continue to use (7o,7jf) = (1)0), and 
to adopt exponential mappings for all radial intervals. 

To verify temporal convergence, we first construct 
an accurate reference solution obtained by evolving the 
perturbed-black-holc initial data to final time tp = 15.0 
with an explicit Dormand Prince 5 (DP5) timestepper 
and timestep At = 0.015625. The spatial domain is de- 
termined by a top spherical harmonic index £ max = 15 
and 1.9 < r < 81.9, and is divided into 8 equally spaced 
concentric shells, each with with N r = 21 radial colloca- 
tion points. Next, for each in a sequence of increasingly 
smaller timcsteps we perform an analogous IMEX evolu- 
tion using the GLoSISDC3 timestepper, which is fourth 
order accurate. One complication involves boundary con- 
ditions: we must ensure that the choices for the explicit 
and IMEX evolutions are consistent. For both we have 
chosen a "frozen" condition, in which the incoming char- 
acteristic is fixed to its initial value, i.e. we freeze U~ b in 



Eq. (|18|) to its initial value. 
We compute the error, 



|A^|| 



max 

a, b 



GLoSISDC3 
b 



DP.5| 
b 



(28) 



and plot it in Figure [3] For intermediate At, we observe 
the predicted fourth-order convergence rate. We remark 
that all timesteps shown in Fig. [5j except the largest, 
correspond to At <C |A| _1 from the standpoint of the 
model problem analyzed in Section IIII Al However, we 
have encountered no stability issues with these short-time 
evolutions. 

We test spatial convergence as follows. Our spatial 
domain, determined by £ max = 15 and 1.9 < r < 41.9, 
is divided into 4 equally spaced concentric shells. For a 
fixed At = 0.0625, we then evolve the perturbed-black- 
holc initial data for different number N r of radial collo- 
cation points in each shell. We compute the root-mean- 
square sum of all constraint violations \[E~ C (see Eq. (53) 
of Ref. [1^] for the precise definition), and plot it in Fig.|6] 
The figure indicates that the solution is dominated by 
spatial error, and exhibits convergence with increased 
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FIG. 5. Temporal convergence test. Error points (circles) 
have been computed using (|28[) in the text. The straight line 
in the plot and its indicated slope have been computed by a 
least squares fit of the third through fifth error points. 
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FIG. 6. Spatial convergence test. This plot depicts histories 
for the constraint energy norm \/Wc described in the text. 



spatial resolution. A plot of the dimensionless constraint 
norm ||C|| denned in Eq. (71) of [25( is qualitatively the 
same. 
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FIG. 7. Stability of various timesteppers when the constraint 
damping terms are treated explicitly or implicitly. Plotted are 
constraint violations vfo- The top two panels show explicit 
treatment of the constraint damping terms. This is stable for 
small timesteps At < 1.024 (top panel) and unstable for large 
timesteps, At > 2.048 (middle panel). The lowest panel shows 
implicit treatment of the constraint damping term, resulting 
in stable evolutions for all timesteps. 



is short enough that the weak instabilities (associated 
with small GLoSISDC3 timesteps) observed in Fig. Q] 
do not arise. Figure [7] shows the constraints for var- 
ious timesteps and three different IMEX timesteppers. 
From the lowest panel, we see that the system is well- 
behaved for all considered timesteps if the constraint- 
damping terms are treated implicitly. The upper two 
panels show that for explicit handling of the constraint 
damping terms, the timestep matters: For small At, the 
simulations behave well, for large At they blow up. This 
is consistent with a Courant limit for the explicit sector 
of the timestepper, arising from the constraint-damping 
term. 



C. Treatment of constraint damping terms 



D. Adaptive timestepping and comparison to 
explicit timestepper 



As described in Sec. Ill A[ the generalized harmonic 
equations ([1} are modified by constraint damping terms 
proportional to 70 in Eq. (|lb[) . These terms cause con- 
straint violations to decay exponentially. Because these 
terms are stiff, they require attention when choosing the 
IMEX splitting, as we now demonstrate. 

We perform runs similar to Fig. [1] but for explicit 
(7jf = 1,7q = 0) and implicit (7^ = 0,7o = 1) con- 
straint damping. The computational domain is the same 
as in Fig. [T]but with Cartesian center (0, 0, 0), N r = 17, 
and L — 9. Our final evolution time for these runs 



In this subsection, we demonstrate adaptive timestep- 
ping in an IMEX evolution by using an adaptive 
timestepper on the perturbed-black-holc initial data from 
Appendix [Cl We evolved this initial data on a set of 16 
concentric spherical shells with Cartesian center (0,0,0) 
and with 1.9 < r < 161.9, N r = 17, and L = 11. A 
gravitational-wave pulse falls into a nonspinning black 
hole of mass M = 1 shortly after t = 0, which causes a 
time-dependent deformation of the hole's horizon. The 
top panel of Fig. [5] shows the minimum and maximum 
values of the intrinsic scalar curvature R of the horizon: 
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FIG. 8. Demonstration of IMEX evolution of a single per- 
turbed black hole using GRrSISDC2 with adaptive timestep- 
ping. Top panel: the minimum and maximum of the horizon's 
dimensionless intrinsic scalar curvature M 2 R, which charac- 
terizes the horizon shape. Bottom panel: The Courant factor 
At/A:E m i n , where At is the size of each timestep and Aa; m j n 
is the minimum spacing between grid-points, for an IMEX 
evolution and for an analogous explicit evolution. Both evo- 
lutions are evolved at the same spatial resolution (with ap- 
proximately 43 3 grid-points). 



As the wave falls into the hole, the horizon shape oscil- 
lates and then relaxes back to the Schwarzschild value 
M 2 R = 1/2, which holds for the curvature of a sphere of 
Schwarzschild radius r = 2M. 

The bottom panel of Fig. [8] plots the step size chosen 
by the adaptive timestepper At/Aa; m i n for an IMEX evo- 
lution and an analogous explicit evolution of the same 
initial data. The explicit timestepper chooses an es- 
sentially constant At, right at its CFL stability limit. 
During the initial perturbation, the IMEX step size de- 
creases to a local minimum; as the hole relaxes to its final 
time-independent configuration, the step size increases, 
eventually reaching an artificially imposed upper limit. 
(This upper limit was chosen to guarantee that the el- 
liptic solver would converge in a reasonable amount of 
wallclock time.) 

During the initial time-dependent perturbation, the 
IMEX evolution is usually able to take significantly larger 
timesteps than the analogous explicit evolution. In 
the explicit evolution, the Courant factor is limited to 
At/ Ax m - m si 3, which is comparable to the minimum of 
the IMEX evolution's Courant factor. 

We remark that the above IMEX simulations exhibit 
some instability: the IMEX run shows slow constraint 
growth, perhaps because we did not impose a constraint- 
preserving boundary condition on the outer boundary. 
However, the analogous explicit evolution exhibits no 
instability, and the IMEX and explicit evolutions' con- 
straint violations are comparable in size when we termi- 
nate the simulations (after time t = 2000M, which is 
long after the spacetime has relaxed to its final, station- 
ary state). 



IV. DISCUSSION 
A. Results obtained in the present work 

In this article, we have further developed IMEX- 
tcchniques applied to hyperbolic systems. Specifically, 
we have moved beyond the model problem of a scalar 
wave [24[ to the study of the full non-linear Einstein's 
equations for single black hole spacetimes. Many results 
of the model problem presented in (2~4| carry over to Ein- 
stein's equations in generalized harmonic form J25|: We 
continue to rewrite the implicit equation in second or- 
der form to utilize an existing elliptic solver [§3. Fur- 
thermore, as in the scalar-field case, we do not impose a 
boundary condition at the excision boundary inside the 
black hole. Uniqueness of the solution of the second order 
implicit equation is enforced, we believe, by the demand 
that the solution be regular across the horizon. 

In contrast to the model problem, the generalized har- 
monic evolution system contains physical constraint^] 
which in explicit simulations are handled with constraint 
damping [J [2^, [2fj[ . We have introduced analogous con- 
straint damping terms in the IMEX formulation, namely 
the terms proportional to Jq in Eqs. ([TU|) and (fTT| . We 
have found that these constraint damping terms are es- 
sential for stability. Treating the constraint damping 
terms explicitly incurs a Courant limit due to their stiff- 
ness, and so we recommend an implicit treatment of these 
terms (7^ = 0;7^ = 70). 

We have focused our investigation on spectral deferred 
correction schemes p7l - l20j , utilizing 3 Gauss-Lobatto and 
2 Gauss-Radau-right quadrature points: GLoSISDC3 
and GRrSISDC2, respectively. These schemes generally 
work well; however, we find a weak instability for small 
timesteps which may be related to the stability region of 
the implicit sector of these IMEX schemes. We also have 
investigated ImexEulcr and third order Additive Rungc 
Kutta (ARK3). While ImexEuler proved robustly sta- 
ble, our simulations with ARK3 showed a linear growing 
instability. The origin of this instability remains an open 
question. 

The most demanding scenario that we have considered 
is a perturbed single black hole that rings down to a qui- 
escent state. We have evolved this configuration with ex- 
plicit and IMEX techniques. The explicit evolution used 
a fifth order Dormand-Prince timestepper with adaptive 
timcstepping; however, because of the necessarily small 
grid-spacing close to the black hole, the explicit simula- 
tion uses an essentially constant timestep at its Courant 
limit, cf. Fig. [5] The IMEX method uses a small timestep 
for the early, dynamic part of the simulation, and then 
chooses increasingly larger timesteps, until it exceeds the 
explicit timestep by about a factor of 200. 



1 These are in addition to the auxiliary constraints arising from 
the reduction to first order form. 
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For very large timesteps, the convergence rate of our 
elliptic solver deteriorates, and overall efficiency drops. 
Therefore, so far we have limited the IMEX timestep to 
« 200 times the explicit timestep. For these timesteps, 
the computational efficiency of the implicit and explicit 
code are approximately similar, for the example shown in 
Fig. [5J We are confident that improved preconditioning 
will accelerate convergence of the implicit solver, allow- 
ing us to utilize yet larger timesteps in IMEX at lower 
computational cost. Besides improved preconditioning, 
several aspects of our future work will increase the effi- 
ciency of the IMEX code: We plan to implement a more 
accurate starting method for the prediction phase of an 
SISDC timestep. We further plan to perform a detailed 
analysis of the required tolerances in the implicit solve 
(in the present work we set tolerances near numerical 
round-off to eliminate spurious instabilities due to insuf- 
ficient accuracy), and we plan to optimize the CH — h code 
implementing Eq. (TIT)) . We expect these steps to signifi- 
cantly increase efficiency of the IMEX code; in contrast, 
the explicit code is already highly optimized. In the next 
subsection, we discuss additional code improvements rel- 
evant to IMEX evolutions of binary black holes. 



B. Prospects for binary black hole evolutions 

Long and accurate binary black hole simulations are 
needed for optimal signal-processing of current and fu- 
ture gravitational wave-detectors |lfl - ll3l ]: this provides 
the motivation for the present work. While the results 
obtained here are very encouraging, additional work will 
be necessary to apply IMEX to black hole binaries. 

First, the formalism must be adopted to the dual- 
frame approach [3lj used in binary black hole simulations 
with SpEC. The corotating coordinates implemented via 
the dual-frame technique are essential for implicit time- 
stepping, because they localize the black holes in the 
computational coordinates. Without corotating coordi- 
nates, the black holes would move across the grid, result- 
ing in rapid time- variability of the solution (on timcscalcs 
M/v, where v denotes the velocity of the black hole with 
mass M). This variability would necessitate a small time- 
step to achieve small time-discretization error. The dual- 
frame technique merely adds a new advection term into 
the evolution equations, therefore, we expect the exten- 
sion to dual-frames to be straightforward. 

Second, the implicit solver must remain efficient de- 
spite the more complicated computational domain. And 
third, good outer boundary conditions will be necessary. 
We expect that the second and third issues can be ad- 
dressed simultaneously with the following ideas: SpEC 
evolves binary black holes on a domain decomposition 
consisting of "inner" spherical shells around each of the 
black holes, which are surrounded by a complicated struc- 
ture of "outer" subdomains (cylinders, distorted blocks 
and spherical shells, the latter of which extend to a large 
outer radius) . The inner spherical shells require the high- 



est resolution and therefore determine the Courant con- 
dition for fully explicit evolutions. 

To simulate binary black holes with IMEX methods, 
we envision a split- by- region approach [32[ , where the in- 
ner spherical shells are treated with the IMEX techniques 
described in this paper and the outer subdomains are 
handled explicitly. The split-by-region approach has two 
important advantages: First, implicit equations will have 
to be solved only on series of concentric shells. This is the 
case considered here, for which SpEC's elliptic solver is al- 
ready reasonably efficient with further possible efficiency 
improvements as discussed in Sec. IIV A"l In contrast, so- 
lution of implicit equations on the entire (rather com- 
plicated) domain-decomposition would likely be less effi- 
cient because of difficulties in preconditioning the inter- 
subdomain boundary conditions. Second, for explicit 
evolutions non-reflecting and constraint-preserving outer 
boundary conditions are available [25], [HI, di[ . Explicit 
treatment of the region near the outer boundary will al- 
low us to reuse these boundary conditions. In contrast, 
similarly sophisticated boundary conditions have not yet 
been investigated in an IMEX setting. 

Because the outer subdomains will be handled explic- 
itly, the split-by-region scheme will still be subject to a 
Courant condition, based on the minimum grid-spacing 
Ax ou ter in the explicitly evolved region. Because the min- 
imum grid-spacing in the outer subdomains is larger than 
the minimum grid-spacing Axi nncl near the black holes, 
the envisioned split-by-region approach should allow for 
timesteps larger by a factor 



R At ee » i. 



(29) 



We shall assume that the cost-per-timestep is propor- 
tional to the number of collocation points, with different 
constants for explicit and IMEX cases: 



Cexplicit — C^outcr + Nij 
ClMEX = CW ou ter + C ' R. 



step -dinner 



(30) 
(31) 



Here, i? s tep is the ratio of the cost of an IMEX-timestep 
to a fully explicit timestep. The simulations presented 
in Sec. Mil give i? s tcp ~ 100, with i? s t ep being somewhat 
larger for very large At and somewhat smaller for small 
At. 

For temporal integration to a fixed final time, the num- 
ber of timesteps for a fully explicit scheme will be propor- 
tional to 1/Axi nner , whereas for the IMEX split-by-region 
scheme, the number of timesteps will be proportional to 
1/Ax outer- Therefore, the IMEX split-by-region scheme 
should require the following fractional amount of CPU 
resources relative to a completely explicit evolution (a 
smaller number indicates advantage for IMEX): 



-Rbbh 



Axw 



IMEX 



1 N n - 



^step dinner 



^explicit RAt 



(32) 
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When -R s tep -dinner > A^outer, this simplifies to 



R 



R 



BBH 



step 



RAt N b 



(33) 



As expected, the question is whether the larger timestcp, 
encoded in i?At, can compensate for the additional cost 
per timestep, encoded in i? s tcp- However, split-by- 
region mitigates the effect of i? s tcp by an extra factor 

Ainner/Atotal- 

To make this discussion concrete, a recent mass-ratio 
q = 6 simulation of non-spinning black holes used N outcl - = 
219222, 7Vinnor= 147288, and R At =34. With these val- 
ues Eq. (|33|) gives -Rbbh = 1-2., i.e. an IMEX evolution 
should be marginally more expensive than a fully ex- 
plicit one. As the mass-ratio is further increased, the 
grid-spacing needed to resolve the smaller black hole de- 
creases proportionally. Therefore, Axi nner will decrease 
proportional to 1/q, and i?At will increase proportional 
to q. The constant of proportionality can be determined 
from i?At = 34 at q = 6, so that i?At ~ 6q. The number 
of grid-points will only modestly change, so we assume 
A^innor ~ Aoutor- Then from Eq. (f3"3"|) we estimate an 
efficiency increase for IMEX of 



-Rbbh 



100 1 
~6q 2 



(34) 



Therefore, with increasing mass-ratio, IMEX will become 
increasingly more efficient than the explicit evolution 
code. 

The additional efficiency gains for IMEX discussed in 
Sec. I IV Al are not taken into account in this estimate. 
Furthermore, a more judicious choice of domain decom- 
position with a more carefully tuned number of colloca- 
tion points in the inner spheres would reduce the ratio 
A'inner/A'totai- Finally, we have not accounted for the fact 
that BBH evolutions require additional CPU resources 
for interpolation. Because interpolation occurs only in 
the outer subdomains, this will reduce i? s tep- 

On the other hand, at this point we do not know how 
accurately the implicit equations must be solved in the 
binary case; if higher accuracy is required to control secu- 
larly accumulating phase-errors, then each implicit solve 
would become more expensive. Furthermore, the binary 
simulations utilize a dual-frame method which will add 
some overhead to the implicit solutions. 

In summary, we believe that IMEX schemes offer the 
promise of faster binary black-hole simulations, but many 
interesting issues (such as those outlined in this section) 
deserve further investigation. 



C. Applicability to other computational techniques 

The results in this paper were obtained for the general- 
ized harmonic formulation of Einstein's equations using 
pseudo-spectral methods. IMEX methods might also be 
implemented for other formulations of the Einstein equa- 
tions, such as the Baumgartc-Shapiro-Shibata-Nakamura 



(BSSN) formulation [3|1 136[ or the recent conformal de- 
compositions of the Z4 formulation J3?J ■ Indeed, for such 
systems specification of the first-order implicit system 
[analogous to Eqs. (flU)) ] corresponding to a single time- 
step is straightforward. However, relative to the analo- 
gous reduction performed for the generalized harmonic 
formulation in this paper, the reduction of such a first- 
order system to a second-order system involving, presum- 
ably, some subset of the system variables would seem to 
be more involved. A second impediment arises from the 
need to use corotating coordinates. In corotating coordi- 
nates, temporal timescales are long, allowing large time- 
steps with sufficiently small time-discretization error (cf. 
Sec. HVB|) . To our knowledge, none of the BSSN/Z4 
codes currently utilize corotating coordinates, although, 
in principle, the dual-frame approach [3l[ could be ap- 
plied in such codes. 

Provided the existence of efficient solvers for the re- 
sulting discretized implicit equations, the IMEX meth- 
ods developed here should also be applicable to other 
spatial discretizations, e.g. finite differences, finite ele- 
ments, or other Galerkin spectral-element approaches. 
The presence of a horizon and the replacement of an 
inner boundary condition by a regularity condition (cf. 
Sec. Ill C| are points demanding particular attention. In 
our approach each component of the apparent horizon 
is covered by a single subdomain. Therefore, in our 
pseudo-spectral treatment the metric in the vicinity of 
the horizon is expanded in terms of a single set of ba- 
sis functions, with regularity of the solution an auto- 
matic consequence. Guaranteed regularity of the solution 
might be lost for either a finite-difference method or an 
unstructured-mesh method, but further studies of these 
possibilities are clearly warranted. 
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Appendix A: Decomposition of T e 

The trace tp ab T eab of the Christoffel symbol T eab of the 
first kind is 

T e = ^da^eb ~ \^ ah d e ^ ab . (Al) 

Writing the time-derivative separately, we reach 

r e =ijj Qb d ip eb + i> kb dki> eb 

- \r h 5° e d^ ab - \r h s k M ab , (A2) 

where is the time t component. Now we insert the 
identities d t ip a b = -MI & + V k $ ka b and d k ip a b = ®kab, 
thereby finding 

r e = -N^ b u eb + ^ b v k $> keb + ^ kb <s> keb 

(A3) 

Finally, we use the identity t a = —Nip 0a to write 

T e = t b U eb - N~ 1 V k t b $ k eb + ^^keb 

- U e 4, ab Ii ab + ^N-^r'te^kab - ^ a Ve®kab- 



(A4) 



Using the last expression, we compute 



Cfl 



- 2^" " N-'vH^r^kab, (A5b) 
and these formulas complete the definitions in Eqs. (1111) . 



Appendix B: Semi-implicit spectral deferred 
corrections 

This appendix describes one of the IMEX timcstepping 
algorithm used for our evolutions, summarizing results 
found in Refs. [l7l - l20j and expressing them in our nota- 
tion. We aim here only to describe the algorithm, and do 
not address stability and convergence issues (which have 
been exhaustively explored in the references). 



1. Collocation approximation of the Picard integral 

We start with the generic ODE initial value problem 
Eq. |4]). Each spectral deferred correction method spec- 
ifies a rule for advancing the vector u n at the present 
timcstep t n (perhaps the initial time to) to a vector u n+ \ 



at the next timestep t n+ i = t n + At. The Picard inte- 
gral form of the initial value problem Eq. (j4|) for starting 
value u„ is 



u(t) = u n + J f(s,u(s))ds. 



(Bl) 



We consider this equation on the interval [i„,t„ + i], 
and show how iterative approximation of (|B1[) yields a 
timcstepping scheme. 

Introduce p collocation nodes which are also time sub- 
steps: 

t( m ) = t„+c m At, < ci < c 2 < • • • < Cp < 1. (B2) 

The c m are either Gauss-Legendre, Gauss-Lobatto, or 
Gauss-Radau nodes relative to the standard interval 
[0,1]. Each of the endpoints, t n and i n +i, may or may 
not be a collocation node. In particular, for the Gauss- 
Legendre case both t n and t n+ i are not in {i(i), • ■ • , t( p )}- 
Define a system vector it( m ) at each collocation point 
£( m ). A solution to the polynomial collocation approxi- 
mation to the Picard integral (|B1[) is a set {«( m ) : m = 
1 , . . . , p} of vectors obeying 



km) 
p 



= U n + At^2 S m qf(t^q), m = l,...,p, 

(B3) 

where ip(t) is the unique degree p — 1 (vector- valued) 
polynomial which interpolates the data 

(i(m),/(i(m),M(m))), m=l,...,p. (B4) 

The elements S mq define the spectral integration matrix. 
The solution {tt( m ) : m = 1, . . . ,p} to Eq. (|B3[) defines 
the approximation 



u, 



tn + l 



if>(t)dt ps u(t n+ i). 



(B5) 



t„ 



We get an approximation to the solution of the colloca- 
tion equations (|B3[) via an iteration described below (that 
is, we get an approximate solution to the approximating 
system of equations). 



2. Iterative solution of the collocation equations 



Our iterative scheme for solving (|B3[) relies on two 
phases: (i) an initial prediction phase which generates 
a provisional solution u® m y and (ii) a correction phase 

which generates successive improvements uk m y k = 
1 , . . . , K as described in detail below. As described by 
Hagstrom and Zhou [2p| . the prediction phase requires 
a starting method, and we use ImexEulcr. For each k 
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the set { u ^ m ^ '■ m — L---,i?} determines an interpo- 
lating polynomial xp k (t), and the numerically computed 
approximation to u(t n+ i) is 



u 



n+l 



U n + 



tn+l 



xl> K (t)dt. 



(B6) 



For the Gauss-Legendre, Gauss- Radau-right, and Gauss- 
Lobatto cases, Hagstrom and Zhou [20| have studied 
the accuracy of these methods. When considered as 
global methods (integration to a fixed time with multiple 
timesteps), they have shown that for sufficiently large K 
the optimal order of attainable accuracy is respectively 
2p, 2p — 1, and 2p ~ 2, that is the same order as for the 
underlying quadrature rule; however, this order is typi- 
cally not obtained for the vectors u¥s at intermediate 
times. 

Typically K = 2p — 1 for Gauss-Legendre, K = 2p — 2 
for Gauss- Radau (left or right), and K = 2p— 3 for Gauss- 
Lobatto cases, where each choice should yield the opti- 
mal order of accuracy. Our presentation of the iteration 
algorithm makes use of the notations 



f(m) 
At Q 
At m 



f{t(m)i u \ m )) 
C X At 

(c m +l - C m ) At, 



(B7) 



1, 



,P~ 1, 



but draws a distinction between two cases (i) Gauss- 
Legendre and Gauss-Radau-right (for these methods t n 
is not a collocation point) and (ii) Gauss-Lobatto and 
Gauss- Radau- left (for these t n is a collocation point). 

To start the prediction phase for Gauss-Legendre and 
Gauss-Radau-right, we first solve 

u il) -At f I (t {1) ,u (1) ) = u n + At Q f E {t n ,u n ) (B8) 

to get For Gauss-Lobatto or Gauss- Radau- left, we 

have v})^ = u(t n ) to start with. We then march forward 
in time by solving in sequence the following equations: 

U °(m+l) — &t m f (i( TO +l)> M ( m +l)) 

= u a {m) +At m f E (t {mhU ° {m) ) (B9) 



for m = 1, . 

defined by the previously constructed w' ml and amounts 
to an ImcxEuler timestep. 

We have used ImexEuler to generate the provisional 



,p — 1. Note that each such equation is 
o 

(m) 



solution {u 



(m) 



= 1, . . . ,p}, and this simple method 
is also the basis of the correction phase. Given { M j>„) : 
m = l,...,p}, a correction sweep yields updated vec- 



tors, {u 



fc+i 
M 



1, . . . ,p\. To understand the eventual 



scheme which produces the updated vectors, first con- 
sider an approximate solution v(t) to the continuum ini- 
tial value problem Eq. (j4}, assuming v(t n ) = u(t n ) = u n , 
along with the residual 



r(t) = u„ 



f(s,v(s))ds - v{t). 



(BIO) 



If the exact solution is u(t) — v(t) + S(t), then the cor- 
rection S(t) obeys 

rift dv 

^L=f(t,v + 5)-f(t,v) + ^, S(t n ) = 0. (Bll) 

We timestep this equation using ImexEuler. For case 
(i), either Gauss-Legendre or Gauss-Radau-right, we first 
solve 

<5(i) - At / / (t(i),'U(i) + 5(i)) = -At / / ^(i),w(i)) +r ( i) 

(B12) 

for S(u , where to reach this equation we have used <5( ) = 
= r(o). For the case (ii) methods we have <5(i) = to 
start with. Subsequently, we solve 

^(m+l) — Atmf 1 {t( m+1 ), ^( m+ i) + 8( m +l)) 
= ^(m) + [f E (t( m ) , V ( m ) + 5( TO ) ) - f E (t (m ) , Vr m } )] 

- Ai m / / (^( m+1 ), U( m+ i)) + rr m+1 ) - r/ m ) 

(B13) 

for m = 1 , . . . , p — 1 . 

To exploit formulas (|B12p — (|B13|) . first make the as- 
signments 

V (m) -> «( m ) 
V( m ) + S {m) u\ m ^ = U k [m) + S k m) 



9=1 



k 

(m)- 



(B14) 



With these assignments in Eq. (|B12[) . we find, upon 
adding w^U to both sides of the equation, 

u*+ 1 -At f I (t (1) ,u k + 1 ) 

^ u (B15) 
= u n -Ai / I (t (1) ,^ 1) ) + At^5i g 4 ) . 

9=1 

Solution of this equation yields . Notice that its 
right-hand side is determined by the known vectors 
{ti^-j : m = 1, . . . ,p}. Next, with the assignments (|B14|) 

in (|B13|) . we find, upon adding u k m+1 -. to both sides, 



U (m+1) &tmf (*(m+l)) u (m+l)) 

= U (m) + (^/ (m) + _ / (m+1) ) 

P 



(B16) 



At ^ Sm+l,qf( q ) — At S rnq f { 



k 

(9)' 



9 =1 g=l 

where we have defined the shorthand 

A/™ = " f E {t (m yu\ m) ). (B17) 
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Sequential solution of this tower of equations yields 
fe+1 for m= 1, ... ,2? — 1. 



u 



(m+l) 



As mentioned, for any method [Gauss-Legendre, 
Gauss- Radau (left or right), Gauss-Lobatto] Eq. (|B6j) de- 
fines a numerical computed approximation to u(t n+ i); 
however, note that for Gauss-Lobatto and Gauss-Radau- 
right, we may also use for this approximation, since 
t( p ) = t n + c p At = t n +i for these cases. 



Appendix C: Perturbed Kerr Initial-Data 

In Sec. IIIII we use initial data representing a nonspin- 
ning Kerr black hole with a superposed gravitational 
wave. Initial data sets are constructed following the 
method of [38j . which is based on the extended confor- 
mal thin sandwich (XCTS) formalism. The Einstein con- 
straint equations read [1^, H(| 



R + K z - KaK 1 



0. 



g lJ K) = 0. 



(CI) 
(C2) 



where Vj is the covariant derivative compatible with the 
spatial metric gij, R = g tJ Rij is the trace of the Ricci 
tensor R^ of g^, and K — g^Kij is the trace of the 
extrinsic curvature Kij of the initial data hypersurfacc. 

The conformal metric gij and conformal factor ip are 
defined by 



g tj = i) 



4^ 



(C3) 



and the time derivative of the conformal metric is denoted 
by 



Uij = (),g,. 



(C4) 



which satisfies Uijg lJ = 0. The conformal lapse is given 
by TV = ip~°N . Applying this conformal decomposition, 
Eqs. (|Cl~|) - (fC2)) can be written as 



1 , ~ 1 
12 



- l^R - ^ b K 2 + -i>- 7 A l3 A lj = 0, (C5) 



( - V, - = 0, (C6) 



and the evolution equation for Kij yields the following 
equation for the lapse: 

V 2 (iV>/> 7 ) - iVV 7 (f + ^ A K 2 + l^A^A^ 

= ^ 5 (d t K-p k d k K). (C7) 

Here (Iffi = V*0» + V j p i - (2/3).^V k f3 k , V, is the 
covariant derivative compatible with g^, R = g 13 Rij 



is the trace of the Ricci tensor Rij of g^, and A %] 



(2N)- 1 ((Lj3)^ 



which is related to Kij by 



Ki 



i>- 10 Aij 



1 



(C8) 



V *J r --vj ■ 

For given g ij: u^, K, and d t K, Eqs. ([U5]). ([06]) . and (|C7|) 
are a coupled set of elliptic equations that can be solved 
for ijj, N, and f3 % . From these solutions, the physical 
initial data g^ and Kij are obtained from ()C3j) and (|C8[) , 
respectively. 

To construct initial data describing a Kerr black hole 
initially in equilibrium, together with an ingoing pulse 
of gravitational waves, we make the following choices for 
the free data, 



9ij 

Uij 

K 
d t K 



4 s 

Ad t hij 

K KS , 
0. 



Ah 



'j • 
1 



~gij~g kl Ad t h 



(C9) 

(CIO) 

(Gil) 
(C12) 



In the above, g^ s and K KS are the spatial metric and 
the trace of the extrinsic curvature in Kcrr-Schild coordi- 
nates, with mass parameter Afxs = 1 an d spin parameter 
o-ks = 0- The pulse of gravitational waves is denoted by 
hij and is chosen to be an ingoing, even parity, m = 2, lin- 



earized quadrupole wave as given by Teukolsky [41|, |42j ) . 
The explicit expression for the spacetime metric of the 
waves in spherical coordinates is 

hijdx L dx 3 = (R\ sin 2 #cos2</>) dr 2 

+ 2i?2 sin 9 cos 9 cos 2cj>rdrd9 

— 2i?2 sin 9 sin 2</>r sin 6drd<f) 

+ [R 3 (l + cos 2 9) cos 2<f> - Ri cos 2</>] r 2 d 2 6 

+ [2 - 2R 3 ) cos 9 sin 20] r 2 sin 9d0d(j) 

+ [-R 3 (1 + cos 2 9) cos 20 + Ri cos 2 9 cos 2<f>] 



'sm*6d 2 



where the radial functions are 



Ri = 3 



Ri 



Rx = - 



F (2) 3F (1) 



^(3) 



2F( 3 ) 



+ 



3F 
6F(i 
9^(2) 



6F 
2LF < 1 

r 4 



(C13) 



(C14) 



(C15) 



+ 



21F 



and the shape of the waves is determined by 
F = F(t + r) = F(x) = e -( x - x ^/ w \ 



d n F(x) 
dx n 



(C16) 

(C17) 
(C18) 



x—t+r 

We choose F to be a Gaussian with width w/M^s = 4 at 
initial radius xo/AIks = 15. The constant A in Eq. (|C9|) 
is the amplitude of the waves. We use the value A = 0.1. 

Equations ([C5]) . (fC6j) . and (jC7j) are solved with the 
pseudospectral elliptic solver described in [3(| ■ 
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